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Abstract 

We present a novel approach for the numerical solution of problems of diffraction by open 
arcs in two dimensional space. Our methodology relies on composition of weighted versions of 
the classical integral operators associated with the Dirichlet and Neumann problems (TE and 
TM polarizations, respectively) together with a generalization to the open-arc case of the well 
known closed-surface Calderon formulae. When used in conjunction with spectrally accurate 
discretization rules and Krylov-subspace linear algebra solvers such as GMRES, the new second- 
kind TE and TM formulations for open arcs produce results of high accuracy in small numbers 
of iterations — for low and high frequencies alike. 

1 Introduction 

Problems of diffraction by infinitely thin open surfaces play central roles in a wide range of problems 
in science and engineering, with important applications to antenna and radar design, electronics, 
optics, etc. Like other wave scattering problems, open surface problems can be treated by means 
of numerical methods that rely on approximation of Maxwell's Equations over volumetric domains 
(on the basis of, e.g., finite-difference or finite-element methods) as well as methods based on 
boundary integral equations. As a result of the singular character of the electromagnetic fields 
in the vicinity of open edges, open-surface scattering configurations present major difficulties for 
both volumetric and boundary integral methods. Boundary integral approaches, which require 
discretization of domains of lower dimensionality than those involved in volumetric methods, can 
generally be treated efficiently, even for high-frequencies, by means of accelerated iterative scattering 
solvers [3, 6, 24]. Unfortunately, integral methods for open surfaces do not give rise, at least in their 
classical formulations, to Fredholm integral operators of the second-kind, and can therefore prove 
to be computationally expensive — as the eigenvalues of the resulting equations accumulate at zero 
and/or infinity and, thus, iterative solution of these equations often requires large numbers of 
iterations. 

This paper presents new Fredholm integral equations of second-kind and associated numerical 
algorithms for problems of scattering by two-dimensional open arcs T (i.e., infinite cylinders of cross- 
section r) under either transverse electric (TE) or transverse magnetic (TM) incident fields. The 
new second-kind Fredholm equations (which result from composition of appropriately modified 
versions and of the classical single-layer and hypersingular integral operators S and N) 
provide, for the first time, a generalization of the classical closed-surface Calderon formulas to the 
open-arc case. In particular, the new formulations possess highly favorable spectral properties: their 
eigenvalues are highly clustered, and they remain bounded away from zero and infinity, even for 



problems of very high frequency. When used in conjunction with spectrally accurate discretization 
rules and Krylov-subspace linear algebra solvers such as GMRES, the new open- arc formulations 
produce results of high accuracy in small numbers of iterations — for low and high frequencies alike. 

The new second-kind formulation for the TM problem is particularly beneficial, as it gives rise 
to order-of-magnitude improvements in computing times over the corresponding weighted hyper- 
singular formulation. Such gains do not occur in the TE case: although the new second-kind TE 
equation requires fewer iterations than the corresponding weighted first kind formulation, the to- 
tal computational cost of the second-kind equation is generally higher in the TE case — since the 
application of the first-kind operator can be significantly less expensive than the application of the 
composite second-kind operator. 

The difficulties that arise as integral equations are used to treat open surface scattering problems 
are of course well known, and many contributions have been devoted to their treatment. Like the 
present work, references [22] and [8] seek to tackle these problems by means of generalizations 
of the classical Calderon relations to the case of open surfaces. The first of these contributions 
establishes that the combination NS can be expressed in the form I + T^, where the kernel 
K(x, y) of the operator has local singularity of at most O ^J_ y ^ - This early result however 
does not take into account the singular edge behavior; the resulting operator Tk is not compact 
(in fact it gives rise to extreme singularities at the edge [16]) and I + Tk is therefore not a second- 
kind operator in any numerically meaningful functional space. When used in conjunction with 
boundary elements that vanish on the edges (by means of well-chosen projections) however, the 
combination NS does give rise to reduction of iteration numbers, as demonstrated in reference [8] 
through numerical examples for low frequency problems. This contribution does not include details 
on accuracy, and it does not utilize integral weights to resolve the solution's edge singularity. A 
related but different method was introduced in [1] which exhibits, once again, low iteration numbers 
for low-frequency problems, but which does not resolve the singular edge behavior and for which no 
accuracy studies have been presented. Finally, high-order integration rules for the single-layer and 
hypersingular operators adapted to open arcs were introduced in a Galerkin framework in [12, 28, 
29]. These methods have thus far only been applied for simple geometries and at low frequencies, 
and limited information is available on the actual convergence properties and performance of their 
computational implementations. 

A second class of methods include those proposed in references [2, 13, 20]. The contribution [2], 
some aspects of which are incorporated in our method, treats the Dirichlet problem for Laplace's 
equation by means of second kind equations. The basis of this approach lies in the observation 
that the cosine basis has the dual positive effect of diagonalizing the logarithmic potential for 
a straight arc and removing the singular edge behavior — so that the inverse of the logarithmic 
potential can be easily computed and used as a preconditioner to produce a second kind operator 
for a general arc. The approach [20], which also uses a cosine basis, treats the Neumann problem for 
the non-zero frequency Helmholtz equation with spectral accuracy by means of first kind equations. 
The contribution [13], finally, treats, just like [2], the Laplace problem by means of second- kind 
equations resulting from inversion of the straight arc logarithmic potential; like [20], further, it 
produces spectral accuracy through use of the cosine transforms. The second-kind integral approach 
developed in [2] and later revisited in [13] seems essentially limited to the specific problem for which 
it was proposed: neither an extension to the Neumann problem nor to the full three dimensional 
problem seem straightforward. And, more importantly, this approach does not lead to adequately 
preconditioned equations for non-zero frequencies: a simple experiment conducted in Section 3.4 
shows that a direct generalization of the algorithm [2] to the Helmholtz problem generally requires 
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significantly more linear algebra iterations than are necessary if the operator S w alone is used. 

The remainder of this paper is organized as follows: after recalling in Section 2 the classical 
boundary integral formulations for the TE and TM open-arc scattering problems, in Sections 3.1 
and 3.2 we present our new weighted operators S w and N w as well as certain periodized counter- 
parts S and N (which are obtained by considering a sinusoidal changes of variables for source and 
observation points). Our main result, the generalized Calderon formula, is presented in Section 3.3. 
Section 3.4 then presents results of numerical evaluation of eigenvalues for a non-trivial open arc 
problem, illustrating the spectral properties of previous open-arc operators as well as the second 
kind operators introduced in this paper. Theoretical considerations concerning the new second- 
kind equations are presented in Section 4, including a succinct but complete proof of the open-arc 
Calderon formulae; a more detailed theoretical discussion, including full mathematical technicali- 
ties, can be found in [16]. The high-order quadrature rules we use for evaluation of the new integral 
operators are described in Section 5. Numerical results, finally, are presented in Section 6 for a wide 
range of frequencies and for various geometries (including a brief study of resonant open cavities) 
demonstrating the uniformly well conditioned character of the integral formulations proposed in 
this paper. 



2 Preliminaries 

Let T denote a smooth open arc in the plane. The TE and TM problems of scattering by the open 
arc T amount to two-dimensional problems for the Helmholtz equation with Dirichlet and Neumann 
boundary conditions on T, respectively: 

Au + k 2 u = outside T, u\ r = f (TE) 

Av + k 2 v = outside T, f^| r = g (TM). ^ 

Here n denotes the normal to T, and / and g are given in terms of the incident electric excitation 
u mc : f = —u mc and g = —du mc /dn on F. 



2.1 Classical open-arc equations for TE and TM problems 

The (unique) solutions u and v of the TE and TM problems above can be expressed, for r T, 
in terms of single- and double-layer potentials of the form u(r) = f T Gk(r, r')/j,(r')d£' and v(r) = 
It 9G dn f ^ v(r')d£'; see e.g. [26,29,30]. Here n r / is a unit vector normal to T at the point r' G T 
(we assume, as we may, that n r / is a smooth function of r' e T), and, calling Hq the first Hankel 
function of order zero, 

GJrA-I iHUk\r-r'\), k>0 

G fc (r,rj- j __^ ln|r _ r/|) k = Q (2) 
and 9G Qnf - = n r / ■ V r 'Gfc(r, r'). Denoting by S and N the single-layer and hypersingular operators 

SM(r)=|G fc (r,r'Mr')rf , r € I\ (3) 
Mr u a 9 f 9G k (r,r') , 

NM(r)= ^y r ^n^ (r) ^ 

d* nm d I dG k (r + zn r ,r>) 
z-+o+ oz J r an r / 



and 



(4) 
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further, the densities [i and v are the unique solutions of the integral equations 

S[m] = /, N[i/] = 5 . (5) 

It is known that the eigenvalues of the integral operators in Equations (5) accumulate at zero and 
infinity, respectively (see e.g. Figure 1). As a result (and as illustrated in Section 6), solution of these 
equations by means of Krylov-subspace iterative solvers such as GMRES generally require large 
numbers of iterations. In addition, as discussed in Section 2.3, the solutions fi and v of equations (5) 
are not smooth at the end-points of T, and, thus, they give rise to low order convergence (and require 
high discretization of the densities for a given accuracy) if standard quadrature methods are used. 
Details in these regards are provided in what follows. 



2.2 Calderon relation for closed and open arcs. 

As is well known [9, 10,21] second-kind formulations for closed surfaces (arcs) result from: (1) The 
classical jump relations associated with singular integral operators, as well as (2) application of the 
Calderon formula 

N c S c = -*+K c , (6) 

which shows that the composition of the closed-surface hypersingular and single-layer operators N c 
and S c equals a compact perturbation K c of the identity. 

In the open-surface context however, jump relations cannot be exploited since, according to 
the boundary conditions (1), the same limits must be reached on both sides of the open surface T. 
Further, use of the combined operator NS does not lead to well-posed equations: as shown in [16], 
the image under NS of the constant function 1 is highly singular, with edge asymptotics NS[l](r) = 

O (^rj)' r G T, where d(r) denotes the Euclidian distance to the edge. The combination NS is 
also problematic in the functional framework set forth in [27, 29, 30] , since the image under S of 
the Sobolev space H2(T) (as defined in [27]) is larger than the domain of definition of H^(T) (as 
defined in [27]) of the operator N. 



2.3 Regularity and singular behavior at the edge 

The singular character of the solutions of equations (5) is well documented [17,27]: \i and v can 
be expressed in the forms 

/i~-^= + r/, v~ X2 Vd + (, (7) 

where d denotes the distance to the edge, xi and X2 are smooth cut-off functions, and where 
the functions r] and £ are somewhat smoother than fi and v. More recently furthermore, it was 
established [17] that expressions such as (7) can be viewed as part of an expansion in powers of 
d 1 / 2 which can be carried out as long as the degree of smoothness of the surface and the right hand 
side of the equations allow it. For example, if the curve itself and the functions / and g in (5) are 
infinitely differentiable, it follows that 

M=^=, v = pVd, (8) 

where a and (3 are infinitely differentiable functions throughout T, up to and including the end- 
points [17]. Thus the singular character of these solutions is fully characterized by the factors d 1 / 2 
and d~ 1 / 2 in equation (8). 
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3 Well-posed Second-Kind Integral Equation Formulations 

3.1 Weighted operators 

In view of the regularity results (8), we introduce a positive integral weight u(r) > with asymptotic 
edge behavior to ~ d 1 / 2 at the edge (by which it is implied that the quotient uj/d 1 / 2 is infinitely 
differentiable up to the edge), and we define the weighted operators 

" ol" 

S w [a] = S - , N w [0]=N[w0], (9) 

-tO J 

so that, in view of the discussion of Section 2.3, for smooth functions / and g, the solutions of the 
equations 

S u [a] = f and N u [fi] = g (10) 

are also smooth. 

Without loss of generality, we choose a smooth parameterization r(t) = (x(t),y(t)) of T defined 
in the interval [—1,1], and we define the weight to canonically as to (r(£)) = Vl -t 2 . As a result, 
the operators S w and N w give rise to the parameter-space operators 

S u [p](t) = fo k (r(£),r(£')) T(t')dt>, (11) 

and 

N^](t) = lim f 1 -j^G k (r(£) + zn r{t) , r(t')) 

z ^o+ dz 7_! <9n r(t ,) ( > ' (12) 

defined on functions ip and ip of the variable t, — 1 < £ < 1; in these equations we have set 
r(£) = |^|. Clearly, for tp(t) = a(r(t)) and Y>(£) = 0(r(t)) we have S w [a](r(t)) = and 
N w [^](r(t)) = JV w [^](t). 

3.2 Periodized spaces and regularity 

Introducing the changes of variables t = cos 9 and t' = cos 6', and defining ng = n r ( cos g), we obtain 
the periodic weighted single-layer and hypersingular operators 

S[j\(0)= G k (r(cos8),r(cos0 / ))~f(8 , )T(cos8 / )d9 / (13) 
Jo 

and 

N[y](9)= lim ^- [ -^-G k (r(cos9) + zn e , r(cos9'))-f{9')T(cos9') sin 2 9'd9'. 

z->0+ OZ Jq OYIqi 

Clearly then, the solutions to the periodic equations 

S[0\ = f, N[4>]=g, (14) 

where f(9) = /(r(cos0)) and g{9) = g(r(cos9)), are related to the solutions of equations (10) by 
the relations 

<p(9) = cp(cos9) , $(0) = ^(cosff). (15) 

In view of (15), the solutions to equations (14) are smooth and periodic, and it is therefore natural 
to study the properties of S and N in the the Sobolev spaces H^(2tt) of 2ir periodic and even 
functions cf. [5, 31]. 
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Figure 1: Eigenvalue distribution for the spiral-shaped arc displayed in Figure 3, with frequency 
j = 100, for the various operators under consideration. Left S. Left-center N. Right-center S^S. 
Right NS. Note the important difference of scale between the four plots. 



Definition 1. Let s£l. The Sobolev space H^(2tt) is defined as the set of 2ir -periodic functions 

oo oo 

of the form v(9) = ^ao + Yl a m cos(m6) for which the s-norm \\v\\g = \ao\ 2 + 2 Yl m2s |om| 2 is 

m=l m=l 

finite. 

Clearly the set {cos(n#) : n £ N} is a basis of the Hilbert space H^(2tt) for all s. 
3.3 Generalized Calderon Formula 

Given the smoothness of the solutions of the equations arising from the weighted periodic operators 
S and N, and in light of equation (6), it is reasonable to consider the composition NS as possible 
basis for solution of open-arc problems. As stated in the following theorem and demonstrated 
below, this combination does indeed give rise to a second-kind integral formulation. The statement 
of the theorem makes use of the concept of bicontinuity: an operator L between two Hilbert spaces 
is said to be bicontinuous if and only if L is invertible and both L and L~ 1 are continuous operators. 

Theorem 1. 

(i) For all s > the operators N and S define bicontinuous mappings S : H^(2ir) — > H^ +1 (2tt) 
and N : H^ +1 (2tt) — > H^(2ir). Thus, the composite operator NS is bicontinuous from H^(2tt) 
into H%(2n) for all s > 0. 

(ii) The generalized Calderon formula 

NS = Jl + K 

holds, where K : H^(2ir) —> H^(2tt) is a compact operator, and where Jq 
is a bicontinuous operator independent of the frequency k. 

(Hi) The set a s of eigenvalues of Jq : H^(2ir) — > H^(2ir) equals the union of the discrete set Aqo = 
{Ao = — ^p, A„ = — | — ^ : n > 0} and a certain open set A s which is bounded away from 
zero and infinity. The sets A s are nested, they form a decreasing sequence, (A s £ A s /,s > s' ) 
and they satisfy P| s>0 A s = { — \}, where A s denotes the closure of A s . 

It follows that the open-arc TE and TM scattering problems (1) can be solved by means of the 
second-kind integral equations 

NS[0\ = N[f] and (17) 



(16) 

: H s e (2n) -> H s e {2ir) 
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NS$] = g, 



(18) 



respectively. The smooth and periodic solutions and ip of these equations are related to the singu- 
lar solutions \x and v of equations (5) via \x (r(cos 6)) = 0(6)/ sin# and v (r(cos0)) = sin 

6(s\f\(6)) ■ 

We thus see that, through introduction of the weight to and use of spaces of even and 27r-periodic 
functions, a picture emerges for the open-surface case that resembles closely the one found for 
closed-surface configurations: the generalized Calderon relation (16) is analogous to the Calderon 
formula (6), and mapping properties in terms of the complete range of Sobolev spaces are recovered 
for S and N, in a close analogy to the framework presented in [14, 21] for the operators N c and S c . 

3.4 Eigenvalue Distributions 

In order to gain additional insights on the character of the various k ^ open-arc operators under 
consideration (namely, S, N, NS as well as a generalization to non-zero frequencies of the operator 
introduced in reference [2]), we consider their corresponding eigenvalue distributions. In Figure 1 
we thus display the eigenvalues associated with these operators, for the spiral-shaped arc displayed 
in Figure 3 and described in Section 6, as they were produced by means of the quadrature rules 
presented in Section 5 and subsequent evaluation of matrix eigenvalues. The frequency was chosen 
to ensure a size to wavelength ratio j = 100, where L is the length of the arc and A the wavelength 
of the incident wave. As expected, the eigenvalues of S tend slowly to zero, the eigenvalues of N are 
large, while the eigenvalues of NS are bounded away from zero and infinity, and they accumulate 
at -\. 

The k > generalization of the equation [2] , whose operator eigenvalues are displayed in the 
center-right portion of Figure 1, is obtained from right-multiplication of the single layer operator 
in equation (10) by the inverse of the flat-arc zero-frequency single-layer operator Sq (defined in 
equation (39) below); the resulting equation is given by 

s&)- i m = f. (19) 

This equation, which can be re-expressed in the form (^1 + (S — <5q)(<So) -1 ^ [0] = / 5 is a second- 
kind Fredholm integral equation: the operator (jS — S^j (Sq)^ 1 : H*(2ir) — > H^(2tt) is compact. 

Unfortunately, the spectrum of the operator in equation (19) is highly unfavorable at high frequen- 
cies, as illustrated in the center-right image in Figure 1. Such poor spectral distributions translate 
into dramatic increases, demonstrated in Table 6, in the number of iterations required to solve (19) 
by means of Krylov subspace solvers as the frequency grows. In fact, a direct comparison with 
Table 3 shows that the second-kind integral equation (19) may require many more iterations at 
non-zero frequencies than the original first-kind equation (14). 

4 Theoretical Considerations 

This section presents a proof of Theorem 1. The argument proceeds by consideration of the flat- 
arc zero- frequency case (Sections 4.1 through 4.3) and subsequent extension to the general case 
(Sections 4.4,4.5 and 4.6). A more detailed discussion, including full mathematical technicalities, 
can be found in [16]. 
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4.1 Flat arc at zero-frequency 

In the particular case where the frequency vanishes (k = 0) and the curve under consideration is 
the flat strip (T = [—1,1]), the operator S reduces to Symm's operator 

S O [0\(O) = -— [ In | cos 6- cos 6'\<p(0')d6', (20) 
27r Jo 

whose well-known diagonal property [5,31] in the cosine basis e n (0) = cosnO, 

( ln_2 n _ q 

SoM = A n e n , \ n = < l_ (21) 

I 2n ' 71 — 1 

plays a central role in our analysis. Note that, in particular, equation (21) establishes the biconti- 
nuity of the operator So from H^(2tt) into H^ +1 (2tt). The corresponding property of bicontinuity 
(from H* +1 (2ir) into H^(2tt)) for the weighted flat-strip, zero- frequency hypersingular operator 

N [tj;]{e) = lim / ^ f ln vV>s0- cosfl') 2 + ^ 2 ) ^(^) sin 2 

results, in turn, from our study of the operator Jo i n Section 4.2. 

An integration by parts argument presented in reference [20] gives rise to the factorization 
N = D S f where 

tMM-^T < 22 > 

and 

f o [0\(9) = ±{<p(O)BmO). (23) 
As a result, for the flat-arc at zero-frequency case, the composite operator NS, is given by 

Jo = N S = D S f S . (24) 

4.2 The operator J 

It is easy to evaluate the action of the operator Jo on the cosine basis: in view of equation (24), 
using (21) and expressing To[e„] as a linear combination of cosines (which results from simple 
trigonometric manipulations) we obtain the relation 

Jo[e n }(d) = \ ft sin rag cosn6> ( 25 ) 

Clearly then, 

j oim = -M - s=!e MW + jf^ ( „)«i„, (26) 

where the operator C is defined by 

C[e n ](9) = \ ° siane 1 = 1 (27) 
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The right hand side of equation (26) may appear to bear a direct connection with the classical 
closed-curve Calderon formula (6). Yet, unlike the operator K c in equation (6), the operator C is 
not compact-see Section 4.3. It is easy, however, to verify that C can also be expressed in the form 



C[0\{9) 



9(ir - 9) 



tt sin ( 



1 



(p{u)du 



7T — ' 



(p{u)du 



(28) 



and is therefore closely related to the Cesaro operator C[rj\(x) = ^ J Q X rj(u)du. Since C is bounded 
(but not compact) from L 2 [0, b] into L 2 [0, b] [4] (where L 2 [0, b] denotes the space of square-integrable 
functions in the interval [0, b]), it follows that C : H^(2tt) — > H*(2tt) is a continuous operator for 
all s > 0, and, therefore, Jo : iJ|(27r) — > H*(2tt) is also a continuous operator. 

Composing the operator Iq = —4S Q ~ 1 CSofo with Jo from the right on one hand, and from the 
left, on the other, and applying the two resulting composite operators to the basis e n , shows that 
Iq is the inverse of Jo, and that Iq is a continuous operator from H^(2tt) into H*(2tt). It follows 
that Jo is a bicontinuous operator. As indicated in Section 4.1, finally, the bicontinuity of Nq 
from iJ| +1 (27r) into H^{2tt) follows directly from equation (24) and the corresponding bicontinuity 
properties of So and Jq. 



4.3 Eigenvalues of J 

Re-expressing (25) in the form 



In 2 

4 

sin(n+l)6> cosng 
4n sin 6 4n 



Jo[e n }(9) 

and making use the well-known expansion [11, 19] 

sin(n + 1)0 



n = 
n>0, 



(29) 



sin# 



= < 



( p 

J2(2-5 0k )cos2k9, n = 2p 

k=0 



we obtain 



Jo[en] 



2 $2 cos(2fc + 1)9, 

k=0 



p-1 



n 



2p + l, 



(30) 



Ane„ - E (1 " ~f)e2k, n = 2p,p>0 

k=0 



1 P_1 
A n e n — 2il 22 e 2k+l, 



(31) 



n 



2p + 1, p > 



k=0 



where the diagonal elements \ n are given by 



In 2 

4 ' 
1 1_ 

4 4n' 



n = 
n > 0. 



(32) 



In view of (29), the operator Jo can be viewed as an infinite upper-triangular matrix. The diagonal 
terms \ n of this matrix are eigenvalues of Jo; the corresponding eigenvectors v n , in turn, can 
be expressed in terms of a finite linear combination of the first n cosine basis functions: v n = 
22k=o c k e k- This establishes that the set Aqo defined in Theorem 1 is contained in the spectrum a s 
for all s > 0. 
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Table 1: Scattering by a spiral-shaped arc of size k = 800: far-field errors. 
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TE case, S 


TE case, NS 


TM case, N 


TM case, NS 


3000 


5.4 x 10~ y 


8.2 x l(r s 


1.5 x IQ-' 6 


2.8 x 10" 4 


3100 


4.5 x 10~ 8 


4.9 x 10~ 8 


1.0 x 10~ 5 


5.5 x 10~ 6 


3350 


8.4 x 10~ 12 


8.5 x 10~ 12 


3.7 x 10" 10 


3.7 x 1Q- 11 



As is known, the diagonal elements of an infinite upper-triangular matrix are a subset of all the 
eigenvalues of the corresponding operator; for instance, the point spectrum of the upper-triangular 

oo 

bounded operator C*[a](n) = (the adjoint of the discrete Cesaro operator C) is given [4, 

k=n 

Th. 2] by the entire disc |A — 1| < 1. A similar situation arises for the operator Jo: searching for 
sequences {f n ,n > 0} such that / = J2 n fn e n satisfies Jo[f] = A/ leads to the relations 



<-i-i>-i = • <33) 

fc=i 

and 



= AA ■ n = - (34) 

k=l 

The spectrum of Jo equals the set of all values A for which the relations (33) and (34) admit 
solutions satisfying ^ \f n n s \ 2 < oo. This set is the union of the discrete set Aoo (see point (iii) in 
Theorem 1), and the open set A s defined by 

A,= {A = (A„A y ),4, + 2< (A -^ ) :|\ g }. (35) 

It follows in particular that the operator C is not compact from H^{2n) into H^(2n): if it were, in 
view of equation (26), the spectrum of Jo would be discrete, in contradiction with (35) 

We thus see that, even though C is not compact (so that, in particular, the decomposition (26) 
does not present Jo as the sum of a multiple of the identity and a compact operator), it follows 
from (32) and (35) that the eigenvalues of Jo are all bounded away from zero and infinity, and that 
they are tightly clustered around —\ — in close analogy to the eigenvalue distribution associated 
with the closed-surface Calderon operator (6). 

4.4 Properties of the operator S 

The study of S in the general case (possibly curved arc, k ^ 0) hinges on the decomposition 

G fc (r(t), r{t')) = A^k, t, t') In \t - t'\ + A 2 (k, t, t'), (36) 
of the kernel Gk, which results from the expression [10, p. 64] 

H^z) = -J (z)ln(z) + G(z), (37) 

7T 

where Jq{z) is the Bessel function and G(z) is analytic. Clearly, the functions Ai and A2 are smooth 
functions of t and t', and the singular behavior of Gk thus resides entirely in a logarithmic term. 



10 



In particular, the continuity of the operator S follows directly from the corresponding continuity 
of the operator Sq. 

To establish the invertibility of the operator S and the continuity of its inverse S" 1 , on the 
other hand, we consider the decomposition 

S = S^I+^y\s-S^, (38) 
where the operator Sq, defined by 

S^] = SoZ [ 7 ], (39) 

differs from So only in the additional multiplicative term ZqItK^) = 7(0)r(cos#). Since the H s e — > 
HI operator (Sq)^ 1 ^ — Sq) is compact (in view of classical embedding results and the increased 
regularity of the kernel of the operator (S — Sq)) and since the operator S is injective (as it follows 
from results established in [27,29,30] for the operator (3)), the bicontinuity of S follows from the 
invertibility of So (and, thus, of Sq) together with an application of the Fredholm theory to the 
term in parenthesis in equation (38). 

4.5 Properties of the operator N 

As is known, for closed curves, the operator N c (the normal derivative of the double-layer potential), 
can be expressed in terms of a composition of a single layer potential and derivatives tangential 
to the curve T c [9, p. 57]. While for open surfaces the conversion of normal derivatives into 
tangential derivatives for the operator N gives rise to highly singular boundary terms (which arise 
from the integration-by-parts calculation that is part of the conversion process), the presence of 
the boundary- vanishing weight oj in our N w operator eliminates all singular terms, and we obtain 
N u [<p](t) = K M (t) + NZ V M (t) where 

KM{t) = k 2 I" G k (r(t), r(0) <ptf) r(t') Vl~t' 2 n t ■ n t , dt' , (40) 

and where 

CMW = ^)J t j\ Gk(r(t),r(t')) A 7!^) dt', (41) 

see [20]. Using the changes of variables t = cos6 and t! = cos#' in equations (40) and (41), it 
follows that 

NM9) = N^](e) + NP v me), (42) 

where 

f'TY 

N 9 \(p\{e) = k 2 / G fc (r(cos^),r(cos^)) r(cos0') 

Jo (43) 

sin 2 9' ne • ne' d0' , 

and where, defining ffi[<p](6) = T ^ 6 ) T o [0\(e), 

Nnm) = v ^ F) (Dosf 5 )m(e). (44 ) 
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Note for future use that equation (44) can be re-expressed in the form 

NV°[<p](6) = N5ffl{0) + -t^-ttMS - SfinmO)- (45) 

where N^} = Z^N Q [7]. 

In view of equation (42), the operator N equals the sum of N 9 (which, like S, maps H^(2tt) 
into H* +l (2tt)) and N pv . But, from (45), we see that N pv can be expressed as the sum of Nq and 
a compact perturbation. Since, like iVo, the operator Nq maps H^ +1 (2tt) into /J|(27r), it follows 
that N pv is a bounded operator from H s e +1 (2-n) into H*(2n), and, thus, N : H s e +1 (2-n) i? e s (2vr) 
is a bounded operator as well. 

The bicontinuity properties of ./V can now be established by an argument similar to the one ap- 
plied to S in the previous section using, this time, the identity N = Nq ( I + (Nq (N — Nq)) . 




Figure 2: Sample diffraction patterns. Left: TE scattering (total field) by an infinitely thin strip 
of size L = 200A for horizontal left-to-right incidence. (The TM total field at this incidence equals 
the incident field.) Right: TM solution for a half-circle of size L = 200A under normal incidence 
from below. Note the trailing shadow in the wake of the strip, and the caustics inside the circular 
reflector. 



4.6 Generalized Calderon formula 

The generalized Calderon formula is obtained by expressing the combination NS in the form 

NS = N^S^ + K = Jq + K (46) 

where the operator K, given by K = N(S — Sq) + (N — Nq)Sq, is compact in view of the classi- 
cal Sobolev embedding theorems and the (increased) smoothness properties of the kernels of the 
operators (S — Sq) and (N — Nq) respectively. In view of the relation Jq = Zq JqZq, Jq and Jq 
have the same spectrum: (X,v) is an eigenvalue-eigenvector pair for Jo if and only if (A, Zq -1 ^]) is 
an eigenvalue-eigenvector pair for Jq. Clearly, Jq has the same mapping properties and regularity 
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as Jo (Jq is bicontinuous from H^(2ir) into H^(2tv)), and therefore, the decomposition (46) shows 
that NS equals the sum of an invertible bicontinuous operator and a compact operator; that is, 
NS is a second-kind Fredholm operator. 



5 High-Order Numerical Methods 

In this section we present spectral quadrature rules for the operators S, N and NS which give rise 
to an efficient and accurate solver for the general open arc diffraction problem (1). 

5.1 Spectral discretization for S 

Use of the nodes |# n = ^^r^ j ; n = 0, - ■ ■ , N — 1, gives rise [23, eq. (5. 8. 7), (5. 8. 8)] to a spectrally- 
convergent cosine representation for smooth, 7r-periodic and even functions ip: 

tf-i ( 2 - Son) N ^ 
(p{9) = ^2 a n cos(n6), where a n = — n ^ <p(9j) cos(n9j). (47) 

n=0 j=0 

Thus, applying equation (21) to each term of expansion (47), we obtain the well-known spectral 
quadrature rule for the logarithmic kernel 

rn N-i 

/ In | cos 9 - cos 9' \<p(9')d9' ~ ^ (p{9 3 )R { p (9) , (48) 

J ° 3=0 

where 

N-l 

R( N \0) = -2 ^ (2 - 5 m )X m cos(m9j) cos(m9). (49) 



m=0 



Following [10, 15, 18] we then devise a high-order integration rule for the operator S, noting first 
from (2) and (36) that 

G fc (r(0), r(0')) = At(k, cos 9, cos 9') In | cos 9 - cos 9'\ + A 2 (k, cos 0, cos 9'), (50) 

where, letting R = |r(cos0) — r(cos6")| we have 

Ai(k, cos 9, cos 9') = - — J (kR), (51) 
2tt 

and 

A 2 (fc,cos0,cos0') = - A Hl{kR) + — J {kR)\n \ cos 9 - cos 9' \. (52) 

4 27T 

In view of (37) and the smoothness of the ratio , — m , the functions A\ and A 2 are even, 

\ ' | COS — COS 0' I ' ± Z ) 

smooth (analytic for analytic arcs) and 27r-periodic functions of 9 and 9' — and, thus, in view 
of (48), the expression 

/ y{9')Ai (jfe, cos 9, cos 0') In | cos 9 - cos 0'|r(cos 9')d9' 
Jo 

N-i (53) 

j=0 
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provides a spectrally accurate quadrature rule. By making use of trapezoidal integration for the 
second term in the right-hand side of (50) we therefore obtain the spectrally accurate quadrature 
approximation of the operator S : 

N-l 

S[<p](0) ~ ^ Yj ( f(0j)r(cose j )(A 1 (k,cose,cose j )Rf\e) + A 2 (k, cos 9, cos 9 j)). (54) 

5.2 Efficient implementation 

The right-hand side of (54) can be evaluated directly for all 9 in the set of quadrature points 
{9 n , n = 0, . . . , N — 1} by means of a matrix-vector multiplication involving the matrix whose 
elements are defined by 

S { n f = ^r(cos^) (A 1 (k,cose n ,cos6 j )Rf ) (9 n ) + A 2 (fc,cos0 n ,cos0j)) . (55) 

A direct evaluation of the matrix on the basis of (49) requires 0(N 3 ) operations; as shown 
in what follows, however, the matrix can be produced at significantly lower computational 
cost. Indeed, expressing the product of cosines in (49) as a sum of cosines of added and subtracted 
angles, the quantities 

N-l 

Rf\0n) = -£(2- 5 m )\ m { cos(^|n - j\) + cos(^(n + j + 1))) (56) 

m=0 

can be expressed in the form 

Rf(9 n ) = rt N H\n -j\)+ rt N \n + j + 1), (57) 
where the vector is given by 

N-l 

R (N) {e) = _J2(2-5 m )\ m cos(™£), £€[0,2N-1]. (58) 

m=0 

Our algorithm evaluates this vector efficiently by means of an FFT, and produces as a result 
the matrix at an overall computational cost of 0(-/V 2 lniV) operations. This fast spectrally- 
accurate algorithm could be further accelerated, if necessary, by means of techniques such as those 
presented in References [3,6,24]. 

5.3 Spectral Discretization for N 

In order to evaluate iV we use (42), the first term of which is a single-layer operator which can be 
evaluated by means of a rule analogous to (54) and a rapidly computable matrix N 9 '^ (similar 
to S( N ^) with elements 

i 2 

N nf N) = ^( cos %) sin2 %K-n e J ^(k, cos 9 n , cos 9 j)Rf\9 n ) + A 2 (k, cos 9 n , cos 9^ . (59) 

To evaluate the second term in (42), in turn, we make use of the decomposition (44), and we 
approximate the quantity Tq [0\ by means of term per term differentiation of the sine expansion of 
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Figure 3: TE-diffraction by a spiral-shaped arc of size L = 200A, for incidence angles of 135° (left) 
and 45° (right) from the positive x-axis. The left figure shows internal reflections that enable the 
field to penetrate to the center of the spiral, giving rise to an interesting array of caustics. The 
relative error e r in both numerical solutions is no larger than 10 -5 . 

TE(5) TE(iVS) 



L 

A 


N 


e r 




Mat. 


It. 


Time It. 


Time 


50 


400 


< 10" 


-5 


< Is 


24 


< Is 8 


< Is 


200 


1600 


< 10" 


-5 


4s 


33 


Is 8 


2s 


800 


6400 


< 10" 


-5 


54s 


45 


18s 8 


15s 



Table 2: Conditioning and computing times for the TE (Dirichlet) problem on the flat strip. 



the function (p{9) sin(0) (which can itself be produced efficiently by means of an FFT). Since Dq is 
essentially the differentiation operator in the x variable, 



1 d 
^9d9 



((^(cos 0)) 



dx 



(*>(*)) 



our solver evaluates the quantity -Do[y>] by invoking classical FFT-based Chebyshev differentiation 
rules [23, p. 195]. 

6 Numerical Results 



The numerical results presented in what follows were obtained by means of a C++ implementation 
of the quadrature rules introduced in Section 5 for numerical evaluation of the operators S and 
./V (and thus, through composition, NS), in conjunction with the iterative linear algebra solver 
GMRES [25]. In all cases the errors reported were evaluated by comparisons with highly- resolved 
numerical solutions. All runs were performed in a single 2.2GHz Intel processor. 
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Table 3: Conditioning and computing times for the TE (Dirichlet) problem on the spiral-shaped 
arc. 
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Table 4: Conditioning and computing times for the TM (Neumann) problem on the strip. 
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Table 5: Conditioning and computing times for the TM (Neumann) problem on the spiral-shaped 
arc. 




Figure 4: A sequence of increasingly thin closed curves converging to the open parabolic scatterer 
x = 1 — 2y 2 and corresponding far-field patterns. From left to right: closed-curve with a=0.9, 
corresponding far field, closed-curve with a=0.99, corresponding far field, parabolic (open) arc 
and corresponding far field. Note the convergence of the far-field patterns as the closed scatterers 
approach the open parabolic scatterer. 
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Spectral convergence. To demonstrate the high-order character of the algorithm described in 
previous sections we consider the problems of TE and TM scattering by the exponential spiral 
x(s) = e s cos(5s), y(s) = e s sin(5s) of size j = 800, where L and A denote the perimeter of 
the curve and the electromagnetic wavelength, respectively. Table 1 demonstrates the spectral 
(exponentially-fast) convergence of the TE and TM numerical solutions produced by means of the 
operators S, N and NS for this problem (cf. equations (14) and (17) and (18)); note from Figure 3 
the manifold caustics and multiple reflections associated with this solution. 

Limit of Closed Curves. In order to obtain an indication of the manner in which an open 
arc problem can be viewed as a limit of closed-curve problems (and, in addition, to provide an 
independent verification of the validity of our solvers) we consider a test case in which the open-arc 
parabolic scatterer x = 1 — 2y 2 is viewed as the limit as a — > 1 of the family of closed curves 
x(s) = (1 — a)coss + acos(2s), y(s) = sin(s). Using the closed-curve Nystrom algorithms [10] 
we evaluate the TE fields scattered by these closed curves at k = 10 for values of a approaching 
a = 1. Figure 4 displays the k = 10 far fields corresponding to a = 0.9 and a = 0.99 side- by-side 
the corresponding far-field pattern for the limiting open parabolic arc as produced by the 5-based 
open-arc solver. Clearly the closed-curve and open-arc solutions are quite close to each other. As 
might be expected, as a approaches 1 an increasingly dense discretization is needed to maintain 
accuracy in the closed-curve solution: for a = 0.9, 256 points where needed to reach a far field error 
of 10~ 4 , while for a = 0.99 as many 1024 points were needed to reach the same accuracy — even 
for the low frequency under consideration. The corresponding open-arc solution, in contrast, was 
produced with 10~ 4 accuracy by means of a much coarser, 64 point discretization. 
Solver performance. The TE (Dirichlet) problem can be solved by means of either the left-hand 
equation in (14) or equation (17), which in what follows are called equations TE(5) and TE(NS), 
respectively. The TM (Neumann) problem, similarly, can be tackled by means of either the right- 
hand equation in (14) or equation (18); we call these equations TM(iV) and TM(NS), respectively. 
Results for TE and TM problems obtained by the various relevant equations for two representative 
geometries, a strip [—1,1] and the exponential spiral mentioned above in this section, are presented 
in Figures 2 and 3 and Tables 2 through 5. In the tables the abbreviation "It." denotes the number 
of iterations required to achieve an e r relative maximum error in the far field (calculated as the 
quotient of the maximum absolute error in the far field by the maximum absolute value of the 
far-field), "Mat." is the time needed to build the S matrix given by equation (55), as well as (when 
required) the corresponding matrix for N 9 which can be constructed in 0{N 2 ) operations from the 
S matrix (see Section 5), and "Time" is the total time required by the solver to find the solution 
once the matrix is stored. 

As can be seen from these tables, the TM equation TM(iV) requires very large number of 
iterations as the frequency grows and, thus, the computing times required by the low-iteration 
second-kind equation TM(NS) are significantly lower than those required by TM(iV). The situation 
is reversed for the TE problem: although, the corresponding second-kind equation TE(iViS') requires 
fewer iterations than TE(5), the total computational cost of the second-kind equation is generally 
higher in this case — since the application of the operator in TE(S') is significantly less expensive 
than the application of operator in TE(iVS). 

Resonant Cavities. We have found that interesting resonant electromagnetic behavior arises 
from diffractive elements constructed as almost-closed open-arcs. As can be seen in Figure 5, for 
example, circular and rocket-shaped cavities with small openings (a few wavelengths in size), can 
give rise to interesting and highly energetic field patterns within the open cavity. The number 
of iterations required for each of these configurations is of course much larger than for simpler 
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Figure 5: Left: TM polarization, for a circular cavity of size j = 200, with a bottom aperture 
of size equal to the wavelength A. Center: TE problem for a rocket-shaped cavity with perimeter 
L = 200A, and a bottom aperture equal to 6A. Right: same as in the center figure, but in the case 
of normal incidence — for which a strong resonance develops in the tail of the rocket. 
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Table 6: Conditioning and times for the solution of the Dirichlet problem using the generalization 
of the method [2, 13]: (5q) -1 as a preconditioner for the spiral-shaped arc. 



geometries, such as the strip. Yet, overall reduction in number of iterations and computing times 
is observed when the equation TM(NS) is used in lieu of TM(iV). For the TE problem, once again 
TE(S') gives rise to faster overall numerics than TE(JVS), although the latter equation still requires 
fewer GMRES iterations. 

Conclusions. We have introduced new second kind equations and associated numerical algorithms 
for solution of TE and TM scattering problems by open arcs. The new open-arc second-kind 
formulations are the first ones in the literature that lead to reduced number of GMRES iterations 
consistently across various geometries and frequency regimes. The new second-kind formulation 
TM(iViS') is highly beneficial for the TM problem, giving rise to order-of-magnitude improvements 
in computing times over the hypersingular formulation TM(iV). Such gains do not occur in the 
TE case: although, the second-kind equation TE(A^S') requires fewer iterations than TE(S'), the 
total computational cost of the second-kind equation is generally higher in this case — since the 
application of the operator in TE(S') is significantly less expensive than the application of operator 
in TE(JVS'). Generalization of these methods enabling efficient second-kind solution of problems of 
scattering by open-surfaces in three dimensions will be presented elsewhere [7]. 
Acknowledgments. The authors gratefully acknowledge support from AFOSR, NSF, JPL and 
the Betty and Gordon Moore Foundation. 
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